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Abstract  (Continued) 


brief  discussion  of  the  motivations  for  developing  such  a global  barotroplc 
model,  describe  In  some  detail  the  numerics  for  the  solution  of  the  model 
equations,  and  present  some  sample  numerical  time-integration  results. 
Several  remarks  concerning  the  numerics  of  this  algorithm  are  also  Included. 
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1.  INTRODUCTION 

The  time  evolution  of  planetary  waves  In  the  troposphere  may  be  described, 
to  a first  approximation,  by  the  barotropic  vorticity  equation  (Rossby  ), 


|f=J<n>>.  <U) 

where  the  absolute  vorticity  n and  the  stream  function  <1/  must  satisfy  the  physical 
constraint 


^ = n - f . (lb) 

Here  f is  the  Coriolis  parameter.  For  a spherical  earth  with  radius  r,  colatitude 
coordinate  o,  longitude  coordinate  X,  and  rotation  rate  fi  about  the  polar  axis,  the 
symbols  In  Eq.  (1)  are  defined  by 


J(rj.  0) 


1 T SlL  Jin.  .M.1 

Sr2slnel308X  3X  30 J * 


(Received  for  publication  9 November  1978) 

1,  Rossby,  C.G.  et  al.  (1939)  Relation  between  variations  In  the  Intensity  of  the 
zonal  circulation  of  the  atmosphere  and  the  displacements  of  the  semi- 
permanent centers  of  action,  J.  Marine  Res.  ,2:38-55. 


1 


[■fe  ain0TW 


V2 3*  - 2 

r sin  0 


sine  ax2 


f = 2fi  cos  0 . 

This  model,  by  virtue  of  its  realism  and  simplicity,  has  been  extensively 

used  in  theoretical  studies  of  the  properties  of  atmospheric  planetary  waves. 

Among  the  large  number  of  such  studies  reported  in  the  meteorological  literature, 

2 

we  may  mention,  for  example,  the  classical  work  of  Fjortoft  on  the  changes  in 

3 

the  spectral  distribution  of  atmospheric  kinetic  energy,  and  the  work  of  Kuo  and 
Lorenz4  on  barotropic  instability.  Indeed,  a model  similar  to  this  one  was  chosen, 
though  with  additional  constraints,  by  Charney  et  al5  to  be  the  first  of  a hierarchy 
of  models  used  to  conduct  numerical  weather  prediction  (NWP)  experiments  for  a 
limited  area  covering  most  of  North  America. 

However,  a solution  in  spherical  coordinates  for  this  model,  simple  as  it  is, 
is  not  easy  to  obtain.  Thus  in  the  case  of  theoretical  studies,  an  infinite  0-plane 
approximation  is  commonly  introduced  to  alleviate  the  geometric  complexity  of 
the  problem.  In  the  case  of  NWP  applications,  this  model  has  been  traditionally 
applied,  using  a polar  stereographic  projection,  only  to  extra-tropical  regions  of 
the  globe.  Solving  Eq.  (1)  for  a limited  region  of  the  globe  requires,  however, 
artificial  lateral  boundary  conditions.  It  is  difficult  and  often  impossible,  to  con- 
struct a set  of  boundary  conditions  which  are  consistent  with  the  differential  equa- 
tions. But  because  the  number  of  arithmetic  computations  needed  for  the  numeri- 
cal solution  of  Eq.  (lb)  on  a sphere  was  prohibitive,  the  dilemma  of  imposing 
extraneous  artificial  lateral  boundary  conditions  to  the  model  equations  was 
accepted  in  the  past  as  a necessary  evil. 

With  the  development  of  computational  mathematics,  it  is  now  possible  to  ob- 
tain an  approximate  solution  to  Eq.  (lb)  on  a sphere  with  relatively  little  computing 
effort  (cf.  Yee6).  It  has  thus  become  less  taxing,  in  terms  of  computer  resources. 


2.  Fjortoft,  R.  (1954)  On  the  changes  in  the  spectral  distribution  of  kinetic  energy 

for  two-dimensional,  non -divergent  flow,  Tellus  5:225-230. 

3.  Kuo,  H. -L.  (1949)  Dynamic  instability  of  two-dimensional  non-dlvergent  flow 

in  a barotropic  atmosphere,  J.  Meteor.  6:105-122. 

4.  Lorenz,  E.N.  (1972)  Barotropic  instability  of  Rossby  wave  motion,  J.  Atmos. 

Set.  29:258-264. 

5.  Charney,  J.  G. , Fjortoft,  R.,  and  von  Neumann,  J.  (1950)  Numerical  integra- 

tion of  the  barotropic  vorticlty  equation,  Tellus  2:2 37-254. 

6.  Yee,  S.  Y.  K.  (1976)  An  efficient  method  for  a finite  difference  solution  of  the 

Poisson  equation  on  the  surface  of  a sphere,  J.  Computational  Phys.  22: 
215-228. 
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to  conduct  controlled  experiments  in  NWP  in  a global  setting.  Among  the  prob- 
lems which  we  plan  to  explore,  using  a global  barotroplc  model,  are  the  following: 
(1)  the  stability  of  barotroplc  flows  on  the  surface  of  a sphere,  (2)  changes  in  the 
spectral  distribution  of  kinetic  energy  in  a barotropic  atmosphere,  (3)  the  effect 
of  non-uniform  spatial  and  temporal  increments  in  a global  model,  and  (4)  the 
effect  of  asstiming  a "free-slip"  boundary  condition  at  the  equator  in  a hemispheric 
model. 


2.  FINITE-DIFFERENCE  APPROXIMATIONS 


If  we  cover  the  surface  of  the  earth  with  a grid  formed  by  the  intersections 
(i,  j)  of  latitudes  and  longitudes,  we  may  approximate  Eq.  (1)  by  the  difference 
analogs 


J i j <n.  0)  , 


(2a) 


(2b) 


where 


j<x.  y>  ■ JA<x,  y)  + ^B(x’  y)  + Jc(x< y)  (2c) 

JA<x.y>  - <xi+l.j  -xi-l(j)<yi,j+l  -yU-l)  -<yt+l.j  -yi-l,j> 

(xi,j+l  -xi.j-l)  - 

JB<x,  y)  “ Iyl,  j+l(xl+l,  j+l  ”xl-l,j+l)  * yl,  j-l(xl+l,  j-l  "xi-l,j-l)1  “ 

fyl+l,  j(xi+l,  j+l  ‘ xl+l,J-l)  ‘ yl-l, j<xi-l,j+l  ' xl-l,  j-l*1  « 

Jc(x.y)  ■ Jxij+i(yi+i,  j+i  -yi-i,j+i>  ‘ xl,J-l(yl+l, J-l  - yi-l.  j-l>1  ‘ 

,xl+l,  j<yi+l,  j+l  "yl+l,j-l)  ‘ xl-l,  j<yi-l,j+l  *yl-l,j-l)1  * 


i 
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♦u*- *1*1-1.  J + 


xu = x(0i' V ■ 


bi  Vi,  j + ci(xi,  j-i  + xi,  j+i*  - (ai  + bi  + 2ci>  \ j 


m.  = 4 AO  AX  sin  0 


sin  0. 


i-f  1^  . 


i-1/2 


AS  sin  0. 


Xi  - j AX  ■ 


AS  = w/I  . 
AX  = 2?r/J  , 


sin  0 


b . 


1+1/2 


2 

AO  sin  0. 


o 2 

AX^  sin  0. 


Here  we  have  divided  a latitude  circle  into  J equal  increments  of  AX  and  a merid- 
ian into  I equal  increments  of  AO.  We  have  also  excluded  the  coordinate  singulari- 
ties at  the  poles  by  constructing  the  grid  in  such  a way  that  the  poles  are  not  grid 
points.  This  is  done  by  placing  the  first  grid-latitude  at  0.  = A0/2.  The  finite- 

7 A 

difference  Jacobian  J . . in  Eq.  (2)  is  due  to  Arakawa,  1 except  that  it  has  been 

J 

modified  here  for  the  spherical  coordinates.  A second  order  five-point  centered- 
difference  formula  has  been  used  for  the  Laplacian 

If  we  now  discretize  the  time  derivative  in  the  left  hand  side  of  Eq.  (2a)  by  a 
centered-difference  "leap-frog"  approximation,  we  have 


(3) 


where  At  is  the  size  of  the  time  increment  and  rj.  . s rj(0{,  X , nAt)  is  the  value  of 

i.  j i j 

n at  grtd-point  (i,  j)  and  time-step  n.  With  these  approximations  in  space  and  time, 
Eq.  (2a)  now  takes  the  form 


n+1  . n-1 

nU  \i 


2 At 
2 

r m . 


J i j (n,  1 1) 


(4) 


Arakawa,  A.  (1966)  Computational  design  for  long-term  numerical  integration 
of  the  equations  of  fluid  motion:  two-dimensional  incompressible  flow,  1, 

J.  Computational  Phya.  1:119-143. 
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For  the  very  first  time-step,  we  use  forward-differencing  in  time, 


(n,  0)  • 


(5) 


By  and  large,  these  are  the  difference  formulae  which  we  shall  use  to  approximate 
Eq.  (1).  We  shall,  however,  experiment  with  the  flexible-increment  approxima- 
tion described  below  for  possible  improvements  in  computational  efficiency. 

One  of  the  more  troublesome  numerical  aspects  in  getting  a solution  for 
Eq.  (2)  is  the  problem  of  computational  instability  due  to  the  convergence  of 
meridians  of  the  spherical  coordinates  at  the  poles.  For  a constant  angular  in- 
crement AX,  the  convergence  of  the  meridians  causes  a decrease,  on  approaching 
the  poles,  of  the  linear  spatial  increment  As  = AXr  sin  0 along  a grid-latitude. 

Thus  in  order  to  maintain  computational  stability,  we  must  use  very  small  At 
values  near  the  polar  regions.  Various  schemes  have  been  proposed  in  the  litera- 
ture to  overcome  this  difficulty. 

In  a previous  study  in  connection  with  the  efficient  finite-difference  solution 

O 

of  Eq.  (2b)  on  a sphere  (Yee  ),  we  adopted  a computational  scheme  in  which  dif- 
ferent sizes  of  angular  increment  are  used  at  different  latitudes  to  approximate  the 
X-derivatives.  We  write  at  latitude  0.,  for  example, 

1 / JLiL 

r sin  I 3X 


. — r sin  0. 
hi  i 


h.  2 AX 


to 


i,  3+hj 


•♦uV 


(6) 


where 


0i,  j+h.  = 0(iA0,  (j  + h.)AX)  . 

Here  h.  is  a positive  integer  which  is  less  than  J and  is  numerically  closest  to  the 
number 

h'  - 1 

ni  ' sin  ei  ‘ 

The  restriction  that  hj  be  integers  is  necessary  here  to  insure  that  the  points 

(i,  j + h.)  coincide  with  the  computational  grid  points.  With  this  formulation,  since 


I 


8.  Yee,  S.  Y.  K.  (1977)  An  Efficient.  Accurate  Numerical  Method  for  the  Solution 

of  a Poisson  Equation  on  a Sphere.  AFGL-TR-77-0246,  Air  Force  Geophysics 
Laboratory,  Hanscom  AFB. 
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t f 


h.  sin  di  ~ 1,  the  linear  increment  As.  = AX  rht  sin  0t  is  approximately  constant 


tor  all  0.. 


This  approach  has  met  with  considerable  success  for  a stable  rapid  solution 
method  for  Eq.  (2b).  We  shall  extend  this  concept  to  the  time-integration  of  Eq. 
(2a).  The  difference  Jacobian  in  Eq.  (2c)  now  takes  the  form 


* ife  * * 

J.  .(x,  y)  s(  JA  (x,  y)  + JB(x,  y)+  Jc<x,  y))/3 


(7) 


where 


JA(x.y>5  (xi+l,  j " xi-l,  j)(yi,  j+h.  'yi,j-hi)  -(yi+l,j  " yi-l,  j)(xi,  j+h(  ' xi,  j-h.* 
* 

JB(x,y)  s (yt>  j+h.<xi+1(  j+h.  " xi-l.j+hi)  ' yi,  j-ht<xi+l.  j-h.  ‘ xi-l,j-hi)l 
Iyi+1,  j<xi+l,  j+hj  ' xi+l,  j-h.'  " yl-l,  j(xi-l,j+h.  ‘ xl-l,J-hl)1  1 

* 

Jc(x.y)  - Ixif  j+h.(yl+1  ( j+h.  - yi-ij+h.)  - xi.  j-hl<yi+i.  j-ht  - yi-i,  j-h"  - 

lxl+l,  j(yl+l,  j+ht  " yi+l,  j-h.'  ‘ xi-l,j(yi-l,j+h.  _ yi-l,j-hl>1  • 

And  Eq.  (4)  becomes 


n+1  n-1  . 2At 

n.  i = 1,  : + 


*n 


U 'U  r2m*  l'j 


J i Ar],ily)  , 


(8) 


where  rcij  = 4 A0  AX  ht  sin  0^  We  reiterate  here  that  with  this  scheme,  time- 
integration  is  to  be  performed  at  every  grid  point  (i,  j).  Only  in  the  difference 
approximations  to  the  X-derivatlves  do  we  invoke  the  idea  of  variable  increments 
in  X. 


3.  RESULTS 

We  have  coded  the  model  Eq.  (2)  in  Fortran  for  test  calculations  in  a CDC 
6600  computer.  In  actual  computations,  we  introduced  the  non-dimensional 
parameters 
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t * tn  , 
i/j*  = 0/(r2n)  , 
n - n/n  , 

so  that  Eq.  (2)  can  now  be  written  as 

=-^~  Jj  j (*!*,  ip*)  . (9a) 

at  mi  *•  J 

i'2  .0*  = n*  - 2 cos  0.  . (9b) 

J * 

This  simple  scaling  reduces  the  amount  of  computation  somewhat  in  our  time- 
integration  procedure.  Grid  resolutions  of  A0  = AX  = 2ir/36  and  At*  - 2»/24 
(At  = 1 hr)  have  been  used  for  the  sample  results  given  in  this  section.  For  initial 
conditions  we  used  variations  of  the  now  classic  Rossby-Haurwitz  waves  given  by 
Phillips,  ® 

ili*  - - Wg  cos  6 + 1*2  cos  6 sinm  6 cos  mX  . (10) 

Here  m is  the  zonal  wave  number  and  w2,  Wg  are  constants  in  units  of  G.  Note 
that  0*  in  Eq.  (10)  increases  with  6 and  is  antisymmetric  about  the  equator.  In 
most  so-called  global  time-integrations,  symmetry  of  the  wind  field  about  the 
equator  is  assumed;  and  time-integrations  are  then  carried  out  only  for  a hemi- 
sphere or  portions  of  a hemisphere.  In  our  time-integrations,  the  following  initial 
fields ‘have  been  adopted: 


(4*  = cos  6 ^ | Wg(m)  sir.m  6 cos  mX(l  - o(m)  - 0.  10(0))  - Wg  cos  0 . (11) 


Here  o(m),  0(0)  are  random  numbers  which  fall  within  the  range  (0,  1).  Note  that 
random  phase-shafts  are  permitted  in  Eq.  (11),  and  the  </>*  field  is  now  no  longer 
antisymmetric  about  the  equator.  For  m2  = rcij  and  a( m)  = 0(0)  = 0,  however,  Eq. 
(11)  reduces  to  Eq.  (10),  which  gives  us  an  unperturbed  initial  field. 

We  present  here  in  Figures  1 through  3 the  stream  function  fields  for  the  cases 
m2  = 1,  2,  and  3 respectively.  In  each  of  these  cases,  nij  is  set  to  equal  1. 


9.  Phillips,  N.A.  (1959)  Numerical  Integration  of  the  primitive  equations  on  the 
hemisphere,  Mon.  Wea.  Rev.  87:333-345. 
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These  figures  are  intended  only  to  give  a feeling  of  the  time  evolution  of  simple 
planetary  waves  in  our  model.  We  have  examined  more  detailed  model  outputs 

and  found  that  for  the  cases  m.  = m„  o(m)  = (3(0)  = 0,  the  phase  velocity  of  the 

1 * 9 

waves  in  our  model  agrees  well  with  that  given  by  Phillips, 


m(m  + 3)w,  - 2Q 

»-  <m  + l)(m  + 2)  ■ <12> 

For  example,  for  the  case  nij  = mg  =1,  Wg  = 0.  06n  , the  waves  indeed  retrograde 
at  a constant  angular  velocity  with  a period  close  to  3.  5 days.  For  the  cases 
presented,  time -integration  was  terminated  arbitrarily  after  30  model  days.  We 
have,  however,  integrated  some  cases  with  unperturbed  initial  fields  (mg  < 3)  for 
up  to  100  model  days  and  detect  no  signs  of  numerical  instability.  Although  not 
needed  to  maintain  computational  stability  for  the  cases  presented  here,  a 4th  order 
linear  filter  designed  by  Shapiro^  has  been  applied  at  each  time-step  to  the  vor- 
ticity  field.  Present  indications  are  that  a Shapiro-type  filter  is  needed  to  main- 
tain numerical  stability  in  long-term  integrations  for  cases  with  mg  > 3.  In  terms 
of  central  processor  time,  a typical  integration  of  2400  time  steps  takes  roughly 
260  sec  in  a CDC  6600.  Thus  both  in  terms  of  stability  and  efficiency,  the  algorithm 
reported  here  seems  to  be  well  suited  for  numerical  experimentation. 


4.  REMARKS 

(1)  Many  studies  have  been  conducted  on  the  numerical  solution  to  barotropic 
models  in  spherical  coordinates.  These  studies  may  be  classified  into  two  broad 
categories:  those  dealing  with  "filtered"  vorticity  equation  models  and  those  deal- 
ing with  "primitive1'  shallow-water  equation  models.  In  the  filtered  equation 

models,  fast  moving  gravity  and  sound  waves  are  "filtered  out"  on  physical  grounds. 

11  12 

Early  numerical  attemps  (for  example,  Silberman,  Baer  and  Platzman  ),  to 
solve  such  model  equations  approximate  the  stream  function  with  a finite  sum  of 
surface  spherical  harmonics,  a natural  choice  for  a basis  function  on  a sphere. 


10.  Shapiro,  R.  (1978)  Personal  communications. 

11.  Silberman,  I.  (1954)  Planetary  waves  in  the  atmosphere,  J.  Meteor.  11:27-34. 

12.  Baer,  F. , and  Platzman,  G.  W.  (1961)  A procedure  for  numerical  integration 

of  the  spectral  vorticity  equation,  J,  Meteor.  18:393-401. 
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This  approach  requires,  however,  a prohibitive  amount  of  arithmetic  computation. 
Thus  only  linearized  flows  or  flows  with  two  or  three  travelling  waves  were 
Investigated. 

The  Initial  success  of  the  barotropic  vortlclty  equation  model  in  operational 

13 

NWP  for  a limited  area  (Cressman  ) apparently  spurred  much  interest  in  the  60's 

14 

to  formulate  global  finite-difference  grids.  The  works  of  Gates  and  Riegel, 

15  16 

Sadourny,  Arakawa  and  Mintz,  0 and  Williamson  are  but  a few  examples  of 
studies  in  this  area. 

With  the  realization  that,  in  terms  of  NWP  skills,  the  quasigeostrophic  model 

has  reached  its  theoretical  plateau  and  that  the  next  higher  order  approximation  in 

filtered  equation  models  involves  mathematical  difficulties,  the  tendency  since  the 

60's  has  been  to  forego  filtering  out  the  fast  waves  and  to  return  to  the  use  of  the 
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"prlmite"  equation  models.  Beginning  with  Phillips,  Kurihara,  Grimmer  and 
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Shaw,  Holloway,  Spelman  and  Manabe,  Williamson  and  Browning,  , Kida, 

22 

and  Umscheid  and  Bannon  among  others  have  contributed  to  the  numerical  study 
of  the  barotropic  primitive  equation  model  on  a sphere.  Unlike  the  simple  scheme 
described  in  Section  2,  the  majority  of  the  work  cited  above  resorted,  however,  to 
complicated  computational  design. 

(2)  With  the  advent  of  the  fast  Fourier  transform,  it  has  become  more  practical 
to  perform  time-integration  of  atmospheric  models  by  the  so-called  pseudo -spectral 


13.  Cressman,  G.  P.  (1958)  Barotropic  divergence  and  very  long  atmospheric 

waves,  Mon.  Wea.  Rev.  86:293-297. 

14.  Gates,  W.L.,  and  Riegel,  C.A.  (1962)  A study  of  numerical  errors  in  the 

integration  of  barotropic  flow  on  a spherical  grid,  J,  Geophys.  Res. 
67:773-784. 

15.  Sadourny,  R.,  Arakawa,  A.,  and  Mintz,  Y.  (1968)  Integration  of  the  non- 

divergent  barotropic  vorticity  equation  with  an  icosahedral-hexagonal  grid 
for  the  sphere,  Mon.  Wea.  Rev.  96:351-356. 

16.  Williamson,  D.  L.  (1968)  Integration  of  the  barotropic  vorticity  equation  on  a 

spherical  geodesic  grid,  Tellus  20:642-653. 

17.  Kurihara,  Y.  (1965)  Numerical  integration  of  the  primitive  equations  on  a 

spherical  grid,  Mon.  Wea.  Rev.  93:399-415. 

18.  Grimmer,  M.,  and  Shaw,  D.  B.  (1967)  Energy-preserving  integrations  of  the 

primitive  equation  on  the  sphere,  Quant.  J.  Roy.  Met.  Soc.  93:337-349. 

19.  Holloway,  J.L. , Spelman,  M.J.,  and  Manabe,  S.  (1973)  Latitude-longitude 

grid  suitable  for  numerical  time  integration  of  a global  atmospheric  model, 
Mon.  Wea.  Rev.  101:69-78. 

20.  Williamson,  D.  L. , and  Browning,  G.  L.  (1973)  Comparison  of  grids  and 

difference  approximations  for  numerical  weather  prediction  over  a sphere, 

J,  Appl.  Meteor.  12:264-274. 

21.  Kida,  H.  (1974)  Tests  of  spherical  grid  systems  for  the  primitive  equations, 

J.  Meteor.  Soc,  Japan  52:1-10. 

22.  Umscheid,  L. , and  Bannon,  P.  R.  (1977)  A comparison  of  three  global  grids 

used  in  numerical  weather  prediction  models,  Mon.  Wea.  Rev.  105:618-635. 
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method  discussed  by  Orszag.  Intense  effort  in  recent  years  in  this  area  of 

research  has  enabled  the  Canadian  Atmospheric  Environment  Service  to  launch 

24 

in  1976  the  first  spectral  model  into  operational  NWP  (Daley  et  al  ).  It  is,  how- 
ever, not  our  purpose  here  to  contrast  the  finite-difference  method  with  the 
pseudo-spectral  method. 

(3)  We  may  look  at  the  differencing  scheme  given  in  Section  2 from  two  dif- 
ferent points  of  view. 

(a)  A non-uniform  angular  increment  AX.  = h.  AX  gives  us  a relatively 
constant  linear  increment  As^  for  all  0^.  Since  h.  increases  toward  the  poles,  we 
have  in  effect,  on  approaching  the  poles,  coarser  and  coarser  resolutions  in  the 
X-derivatives.  The  time-integration  performed  at  every  grid  point  is  really 
nothing  more  than  just  a means  by  which  we  avoid  altogether  the  problem  of  inter- 
polation usually  associated  with  a grid  of  mixed  resolution. 

(b)  The  use  of  a variable  AX^  for  the  X-derivatives  may  also  be  interpreted 
as  the  use  of  a variable  Atj  in  Eq.  (8),  for  we  may  write  Eq.  (8)  as 


where  At.  = At/h..  Since  h^  ~ 1/sin  0.,  AC  decreases  toward  the  poles.  Thus  the 

poleward  increase  of  AXj  in  the  evaluation  of  X-derivatives  is  equivalent  to  keeping 

a uniform  AX  but  reducing  the  size  of  At.  on  approaching  the  poles.  This  is  some- 

l 1ft 

what  similar  to  the  method  given  by  Grimmer  and  Shaw. 

(4)  For  our  flexible  differencing  scheme,  the  truncation  error  in  X-deriva- 
tives increases  toward  the  poles.  This  may  turn  out  to  be  a problem  for  certain 
flow  patterns  as  far  as  accuracy  is  concerned.  This  potential  problem  area 
deserves  and  will  receive  our  closest  attention. 


23.  Orszag,  S.  A.  (1970)  Transform  method  for  calculation  of  vector  coupled  sums: 
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24.  Daley,  R.,  Girard,  C.,  and  Henderson,  J.  (1976)  Short-term  forecasting  with 
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